Modeling elastic instabilities in nematic elastomers 
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Liquid crystal elastomers are cross-linked polymer networks covalently bonded with liquid crystal 
mesogens. In the nematic phase, due to strong coupling between mechanical strain and orient at ional 
order, these materials display strain-induced instabilities associated with formation and evolution 
of orientational domains. Using a 3-d finite element elastodynamics simulation, we investigate one 
such instability, the onset of stripe formation in a monodomain film stretched along an axis per- 
pendicular to the nematic director. In our simulation we observe the formation of striped domains 
with alternating director rotation. This model allows us to explore the fundamental physics gov- 
erning dynamic mechanical response of nematic elastomers and also provides a potentially useful 
computational tool for engineering device applications. 

PACS numbers: 61.30.Vx, 62.20.D-, 64.70. mf, 81.40.Jj, 83.80.Va, 83.80.Xz 
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Liquid crystal elastomers (LCE) exhibit some of the 
elastic properties of rubber along with the orientational 
order properties of liquid crystals, displaying a variety of 
nematic and smectic phases. They are composed of liq- 
uid crystal mesogens covalently bonded to a cross-linked 
polymer backbone [l|, [2[ . These materials display strong 
coupling between orientational order of the mesogens and 
mechanical deformation of the polymer network. For in- 
stance in a nematic LCE, any change in the magnitude of 
the nematic order parameter induces shape change, e.g. 
the isotropic-nematic phase transition induces strains of 
up to several hundred percent [3[. 

Conversely, applied strain can also drive changes in ori- 
entational order, producing the fascinating phenomenon 
of semisoft elasticity [4]. In a classic experiment, Kundler 
and Finkelmann [5[ measured the mechanical response of 
a monodomain nematic LCE thin film stretched along 
an axis perpendicular to the nematic director. They ob- 
served a semi-soft elastic response with a pronounced 
plateau in the stress-strain curve arising at a threshold 
stress. Accompanying this instability they observed the 
formation of striped orientational domains with alternat- 
ing sense of director rotation, and a stripe width of 15 
fim. They repeated the experiment with samples cut at 
different orientations to the director axis, and found that 
the instability was absent when the angle between the 
initial director and the stretch axis was less than 70° ; 
in this geometry, instead of forming stripes, the director 
rotates smoothly as a single domain. 

DeSimone et al 0] carried out numerical simulation 
studies of the stripe instability using a two-dimensional 
finite element elastostatic method. Each area element 
in the system was considered as a composite of domains 
with different orientations. This simulation model was 
the first to reproduce successfully the soft elastic response 
of nematic elastomers, but did not attempt to resolve the 
resulting microstructural evolution. 

Here we explore this elastic instability in more detail 



by simultaneously modeling the sample's mechanical re- 
sponse and the associated microstructural evolution as 
a function of strain. We use a Hamiltonian-based 3-d 
finite element elastodynamics model with terms that ex- 
plicitly couple strain and nematic order. By resolving the 
finite element mesh down to the micron scale, we resolve 
the formation of orientational domains, and because the 
model is dynamic rather than static in character, we can 
examine the effects of strain rate. We use the simula- 
tion to explore the dependence of mechanical response 
on deformation geometry. 

We model this instability in a thin film of nematic 
elastomer which has been cross-linked in the nematic 
phase. Using public domain meshing software |7| we 
discretize the volume of the sample into approximately 
78, 000 tetrahedral elements. For each volume element 
we assign a local variable n which defines the nematic 
director, and Qij = ^S(3niTij — Sij) which is the associ- 
ated symmetric and traceless nematic order tensor. The 
initial state is taken to be a monodomain with n = n Q 
in every element; this configuration is defined as the sys- 
tem's stress-free reference state. 

There are many approaches to finite element simu- 
lation of the dynamics of elastic media [8]; we make 
use of an elegant Hamiltonian approach developed by 
Broughton et al (§, Qij, generalizing it to three dimen- 
sions and the case of large rotations. We write the Hamil- 
tonian of an isotropic elastic solid as: 



lastic 
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(1) 



Here the first term represents elastic strain energy, with 
p summing over volume elements. V p is the volume 
of element p in the reference state. For an isotropic 
material the components of the elastic stiffness tensor 
Cijki are determined from only two material parame- 
ters, namely the shear and bulk moduli [ll|. As an 
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approximation, Broughton et al developed this formu- 
lation using the linear strain tensor, but we instead use 
the rotationally invariant Green-Lagrange strain tensor 



u k,i u k,j), where u is the displace- 



ment field. We note that using the linearized strain ten- 
sor would make the Hamiltonian unphysical, as rotation 
of the sample would appear to cost energy. The sec- 
ond term represents kinetic energy in the lumped mass 
approximation whereby the mass of each element 
is equally distributed among its vertices, which are the 
nodes of the mesh. Here i sums over all nodes, mi is the 
effective mass and vi the velocity of node i. 

To account for the additional energy cost associated 
with the presence of a director field, we add to the po- 
tential energy, 



H nematlc = Y,{-<Mi ~ + / Wi - QW 2 ] 
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The first term describes coupling between the strain and 
order parameter tensors using a form proposed by De- 
Gennes [12[. Here Q?- defines the nematic order in the 
element's reference state. The prefactor a controls the 
strength of this coupling, and DeGennes [13] argued that 
it is of the same order of magnitude as the shear mod- 
ulus \i. Variables Qr/, Q°p and Sij are all defined in 
the body frame, i.e. they are invariant under rotations 
in the target frame. See [13] for the relation between 
Qij in the body and lab frames. The second term de- 
scribes "cross-link memory," that is, the tendency of the 
nematic director to prefer its orientation at crosslinking. 
Thus there is an energy cost to rotate the director away 
from its reference state, with coupling strength (3. The 
third term is an energy penalty for spatial variations of 
the nematic director, similar to a Frank free energy in 
the single elastic constant approximation. The summa- 
tion is carried only over nearest neighbour elements in 
the mesh, as the typical domain size is of the order of the 
nematic correlation length [14]. 



The strain tensor within each tetrahedral element 
is calculated in two steps. We calculate the displacement 
u of each node from the reference state, then perform 
a linear interpolation of the displacement field within 
the volume element in the reference state. The result- 
ing interpolation coefficients represent the derivatives Uij 
needed to calculate the components of the strain tensor. 
Details can be found in any introductory text on finite 
element methods, e.g. [T^j. At this level of approxima- 
tion, the strain is piecewise constant within each volume 
element. The effective force on each node is calculated 
as the negative derivative of the potential energy with 
respect to node displacement. To include internal dissi- 
pation in the system, we add an additional force term 





FIG. 1: (Color online) Simulation: stretching a nematic elas- 
tomer film at an angle of 90° to the director. Initially a 
monodomain, the director field evolves to form a striped mi- 
crostructure. 



which depends on the local velocity gradient [l6j |. This 
dissipation is isotropic in character and does not depend 
on the orientation of the director field. 

To evolve the system forward in time, we assume the 
director is in quasistatic equilibrium with the strain; that 
is, the time scale for director relaxation is much faster 
than that for strain evolution as observed by Urayama 
[l 71 ] . The first part of each step is elastodynamics: hold- 
ing Qij in each element constant, the equations of mo- 
tion / = ma for all node positions and velocities are 
integrated forward in time using the Velocity Verlet al- 



gorithm [18|], with a time step of 10 sec. In the second 



part of each step, we relax the nematic director in each 
element to instantaneously minimize the element's poten- 
tial energy. Because the director rotates from a higher 
energy state to a lower energy state without picking up 
conjugate momentum, this is a source of anisotropic dissi- 
pation. Thus in our model, as in real nematic elastomers, 
strains that rotate the director cause more energy dissi- 
pation than those applied parallel to the director [19[ . 

We simulate uniaxial stretching in an initially mon- 
odomain nematic elastomer film of size 1.5 mm x 0.5 
mm with a thickness of 50 /im, with shear modulus 
fj, = 5.7 x 10 5 Pa, bulk modulus B r = 2.8 x 10 7 Pa, 
and parameters a = /i, (3 = 0.3/i, and 7 = 10 -7 . We 
first consider the case where the director is initially ori- 
ented along the y axis, transverse to the direction of ap- 
plied strain. The sample is clamped on two sides and 
the clamped regions are constrained to move apart lat- 
erally at a constant speed of 1 mm/sec. The resulting 
microstructural evolution is shown in FigQ3 Here color 
represents Jones matrix imaging of the director field as 
viewed through crossed polarizers parallel to the x and y 
directions; blue corresponds to a director parallel to the 
polarizer or analyzer, and red corresponds to a director 
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FIG. 2: (Color online) Engineering stress (circles) and direc- 
tor rotation (squares) vs applied strain, for the system shown 
in Fig. 90degrees. Onset of director rotation and the stress- 
strain plateau both occur at the same strain. 



at a 45° angle to either. While the simulated sample 
is three-dimensional, the film's microstructure does not 
vary significantly through the thickness and can thus be 
visualized in 2-d. 

At a strain of 8.5%, the director field in the sample be- 
comes unstable and orientational domains form, nucleat- 
ing first from the free edges of the film. By 9% strain, the 
whole film is occupied by striped orientational domains 
with alternating sense of director rotation. The stripes 
are not uniform in width, being slightly larger near the 
free edges. Near the center of the sample, each individual 
stripe has a width of about 25 /xm, which is of the same 
order of magnitude as that observed in experiment [5[. 
This value is in reasonable agreement with the theoreti- 
cal estimate by Warner and Terentjev [1] who predicted 
a stripe width of h ~ y/^L/yl — 1 / Af ; where £ is the 
nematic penetration length, L is the sample width, and 
Ai is the strain threshold of the instability. The stripes 
coarsen as the elongation increases. Eventually this mi- 
crostructure evolves into a more disordered state with 
stripes at multiple orientations. By reaching 35% strain, 
the stripes have vanished and the film is again in a mon- 
odomain state with the director oriented with the direc- 
tion of strain. Only the regions near the clamped edges 
do not fully realign, in agreement with experimental ob- 
servations [5[ and with the simulation studies of DeSi- 
mone [6| . We will explore the dependence of stripe width 
on aspect ratio and other parameters in future work. 

The resulting stress-strain response is semi-soft [HI in 
character, as shown in FigJ21 The initial elastic response 
is linear, followed by an extended plateau running from 
about 8.5% to over 30% strain, after which there is a 
second linear regime. We also measure the average direc- 
tor rotation < sin 2 (aV) > and observe that the thresholds 
for both the stress-strain plateau and the rotation of the 
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FIG. 3: (Color online) Engineering stress (circles) and direc- 
tor rotation (squares) vs applied strain, applied at an angle 
of 60° from the nematic director. 





FIG. 4: (Color online) Simulation: stretching a nematic elas- 
tomer film at an angle of 60° to the director. Initially a mon- 
odomain, the director field rotates smoothly without sharp 
gradients in orientation. 



nematic director occur at the same strain. This finding 
demonstrates, in agreement with theory [H, HH, that the 
reorientation of the system's internal degree of freedom- 
namely the nematic director-reduces the energy cost of 
the deformation. 

We also performed simulations for monodomain ne- 
matic elastomer films with the initial director orientation 
at different angles to the pulling direction. In Figj3] we 
plot the film's stress-strain response when strain is ap- 
plied at an angle of 60° from the nematic director, which 
shows no plateau, and likewise director rotation shows no 
threshold behavior. As shown in FigJH the director ro- 
tates to align with the strain direction without forming 
stripes. We performed additional simulations with the 
director at angles of 70° and 80° to the pulling direction 
and again found no stripe formation and no plateau in 
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FIG. 5: (Color online) Dependence of the stress-strain re- 
sponse on strain rate. 



will also examine the role of initial microstructure and 
thermomechanical history in determining mechanical re- 
sponse. Using the same finite element approach, we can 
also test the predictions of other proposed constitutive 
models, and model geometries of interest for potential 
applications. Through this approach we hope to bridge 
the divide between fundamental theory of these fascinat- 
ing materials and engineering design of devices. 

We thank Profs. T. Lubensky, M. Warner, E. Terent- 
jev and A. Desimone for fruitful discussions. This work 
is supported by NSF DMR-0605889, daytaOhio, and the 
Ohio Supercomputer Center. FY also acknowledges the 
Institute for Complex Adaptive Matter - Branches Cost 
Sharing Fund for a postdoctoral fellowship. 




FIG. 6: Simulation: A nematic elastomer disk is stretched 
radially. The director field smoothly transforms from a ho- 
mogeneous monodomain to a radial configuration. 



the stress-strain response. 

We also tried varying the applied strain rate. Fig|5] 
compares the stress-strain response for samples strained 
at 1 mm/sec and 5 mm/sec. The higher strain rate pro- 
duces a significant stress overshoot, and stripe formation 
occurs at a strain of 15%. This finding suggests that the 
threshold strain for the instability depends in a signifi- 
cant way on strain rate. 

Next we simulated a monodomain nematic elastomer 
film under isotropic strain. A circular sample of diame- 
ter 1 cm and thickness 100 /im, with the director initially 
along the y axis, was stretched radially in all directions 
by pulling the edge outward at constant speed. Figj6] 
shows the film at different stages of its extension, demon- 
strating that the director field smoothly changes from a 
monodomain to a radial configuration. With a careful 
choice of the sample's thickness, this deformed circular 
sheet of nematic elastomer could be used as a tunable 
spatial polarization converter as described in [20[. 

The simulations presented here were performed at far 
higher strain rates, e.g. 50% per second, than those used 
in typical experiments [5|, [2l| where the material is al- 
lowed to relax for minutes or hours between strain in- 
crements. In future work we plan to apply our model 
to examine deformation of nematic elastomers at slower 
strain rates and as a function of sample geometry. We 
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